rm(list = ls())
require(miceadds)
require(stats4)
require(lmtest)
require(bbmle)
require(msm)
require(nlWaldTest)
require(sandwich)
require(stargazer)
require(mfx)
require("multiwayvcov")
require(car)
require(alr3)
 
#data2old=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/infoexp2_data/responses.csv", header=TRUE)
#data13=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_july25_type13clean.csv", header=TRUE)
#data15=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_August12_type15(clean).csv", header=TRUE)

data2old=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/infoexp2_data/responses.csv", header=TRUE)
data13=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_july25_type13clean.csv", header=TRUE)
data15=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_August12_type15(clean).csv", header=TRUE)


dataframe2old=data.frame(data2old)
dataframe13=data.frame(data13)
dataframe15=data.frame(data15)


priors=c(0.5,0.6,0.75,0.85)

expframe=rbind(dataframe13,dataframe15,dataframe2old)
expframe=subset(expframe,uid>200)

#note that player 234 did not finish and so has no real data
expframe=subset(expframe,uid!=234)


#definitions
expframe$priorA=(expframe$qid==1)*0.5+(expframe$qid==5)*0.6+(expframe$qid==6)*0.75+(expframe$qid==7)*0.85
expframe$priorB=1-expframe$priorA
expframe$Aresponse=(expframe$response==6)
expframe$Bresponse=(expframe$response==7)
expframe$stateA=(expframe$state==4)
expframe$stateB=(expframe$state==10)
expframe$correct=expframe$stateA*expframe$Aresponse+expframe$stateB*expframe$Bresponse
uidlist=unique(expframe$uid)

#find N
#**
length(uidlist)
#**

#Begin constructing first table
pAgiven2=vector('numeric')
restrictionpAgiven1=vector('numeric')
pAgiven1=vector('numeric')
prob=vector('numeric')
sterr=vector('numeric')
tstat=vector('numeric')


#Base Regression
modelCP=lm(Aresponse~factor(stateA+priorA)+0,data=expframe)
summary(modelCP)



#Adjust COV
vcovmodCP=vcovHC(modelCP, type="HC3")
vcovmodCPcluster=cluster.vcov(modelCP, expframe$uid)

modelcpsterr=vector("numeric")
#Model CP sterrors
for(i in 1:8){
modelcpsterr[i]=sqrt(vcovmodCPcluster[i,i])
}

#Experiment 3 Regression logit
logitCP=glm.cluster(data=expframe, Aresponse~factor(stateA+priorA)+0, cluster="uid", weights=NULL, subset=NULL,  family="binomial" )



pAgiven2[1]=modelCP$coefficients[1]
pAgiven2[2]=modelCP$coefficients[2]
pAgiven2[3]=modelCP$coefficients[3]
pAgiven2[4]=modelCP$coefficients[4]

pAgiven1[1]=modelCP$coefficients[5]
pAgiven1[2]=modelCP$coefficients[6]
pAgiven1[3]=modelCP$coefficients[7]
pAgiven1[4]=modelCP$coefficients[8]

restrictionpAgiven1[1]=((2*priors[1]-1)/priors[1])+pAgiven2[1]*((1-priors[1])/priors[1])
restrictionpAgiven1[2]=((2*priors[2]-1)/priors[2])+pAgiven2[2]*((1-priors[2])/priors[2])
restrictionpAgiven1[3]=((2*priors[3]-1)/priors[3])+pAgiven2[3]*((1-priors[3])/priors[3])
restrictionpAgiven1[4]=((2*priors[4]-1)/priors[4])+pAgiven2[4]*((1-priors[4])/priors[4])

sterr[1]=as.numeric(deltamethod(~x5-((2*0.5-1)/0.5)+x1*((1-0.5)/0.5),modelCP$coef,vcovmodCPcluster))
tstat[1]=(pAgiven1[1]-restrictionpAgiven1[1])/sterr[1]
prob[1]=1-pnorm(tstat[1])

sterr[2]=as.numeric(deltamethod(~x6-((2*0.6-1)/0.6)+x2*((1-0.6)/0.6),modelCP$coef,vcovmodCPcluster))
tstat[2]=(pAgiven1[2]-restrictionpAgiven1[2])/sterr[2]
prob[2]=1-pnorm(tstat[2])

sterr[3]=as.numeric(deltamethod(~x7-((2*0.75-1)/0.75)+x3*((1-0.75)/0.75),modelCP$coef,vcovmodCPcluster))
tstat[3]=(pAgiven1[3]-restrictionpAgiven1[3])/sterr[3]
prob[3]=1-pnorm(tstat[3])

sterr[4]=as.numeric(deltamethod(~x8-((2*0.85-1)/0.85)+x4*((1-0.85)/0.85),modelCP$coef,vcovmodCPcluster))
tstat[4]=(pAgiven1[4]-restrictionpAgiven1[4])/sterr[4]
prob[4]=1-pnorm(tstat[4])

NIASframe=data.frame(priorA=priors,pagiven2=pAgiven2,restriction=restrictionpAgiven1,pagiven1=pAgiven1,prob=prob)

#*****************************
stargazer(NIASframe,rownames=FALSE,summary=FALSE)
#*****************************

